setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
sink("log_replication.R")

rm(list=ls(all=TRUE))
library(lm.br)
library(tidyverse)
library(MCMCpack)
library(stargazer)
library(ggplot2)
library(strucchange)
library(mcp)
library(ecp)
library(readstata13)
library(dplyr)
library(ranger)
library(xtable)
library(pROC)

unilateral_preds = read.csv("unilateral_preds.csv")

## prepare data
unilateral_preds$year = as.numeric(as.character(unilateral_preds$year))
unilateral_preds$doctype[unilateral_preds$doctype == "eo"] = "EO"
unilateral_preds$doctype[unilateral_preds$doctype == "pr"] = "PR"
unilateral_preds$doctype = car::recode(unilateral_preds$doctype, "c('04','21','53','05','08','12','17','20','37','52','56') = 'MM';c('EO', '03', '06', '22','33','41') = 'EO'; c('29', 'PR', '35') = 'PR'")
unilateral_preds = unilateral_preds[unilateral_preds$doctype %in% c("EO", "MM", "PR"),]
unilateral_preds$doctype = as.factor(as.character(unilateral_preds$doctype))
unique(unilateral_preds$doctype)
unilateral_preds = unilateral_preds %>% filter(year >= 1877 & year < 2021)

## create issues labels
unilateral_preds$issue[unilateral_preds$policyarea_1==1] = "Economic"
unilateral_preds$issue[unilateral_preds$policyarea_1==2] = "Civil rights"
unilateral_preds$issue[unilateral_preds$policyarea_1==3] = "Health"
unilateral_preds$issue[unilateral_preds$policyarea_1==4] = "Agriculture"
unilateral_preds$issue[unilateral_preds$policyarea_1==5] = "Labor"
unilateral_preds$issue[unilateral_preds$policyarea_1==6] = "Education"
unilateral_preds$issue[unilateral_preds$policyarea_1==7] = "Environment"
unilateral_preds$issue[unilateral_preds$policyarea_1==8] = "Energy"
unilateral_preds$issue[unilateral_preds$policyarea_1==9] = "Immigration"
unilateral_preds$issue[unilateral_preds$policyarea_1==10] = "Transportation"
unilateral_preds$issue[unilateral_preds$policyarea_1==12] = "Law/Crime"
unilateral_preds$issue[unilateral_preds$policyarea_1==13] = "Social welfare"
unilateral_preds$issue[unilateral_preds$policyarea_1==14] = "Housing"
unilateral_preds$issue[unilateral_preds$policyarea_1==15] = "Domestic commerce"
unilateral_preds$issue[unilateral_preds$policyarea_1==16] = "Defense"
unilateral_preds$issue[unilateral_preds$policyarea_1==17] = "Technology"
unilateral_preds$issue[unilateral_preds$policyarea_1==18] = "Trade"
unilateral_preds$issue[unilateral_preds$policyarea_1==19] = "International affairs"
unilateral_preds$issue[unilateral_preds$policyarea_1==20] = "Government operations"
unilateral_preds$issue[unilateral_preds$policyarea_1==21] = "Public lands"


## limit to significant directives only
significant_preds = unilateral_preds %>%
	filter(Significance_rf >= .355)

## create annual measures
unilateral_annual_type = unilateral_preds %>%
			 		group_by(year, doctype) %>% summarize(numdocs = n())


## create annual measures, limited to significant directives
significant_sum = unilateral_preds %>% 
				mutate(is.sig = Significance_rf >= 0.355) %>%
				group_by(year) %>% 
					summarize(numsig = sum(is.sig),
							numsigEO = sum(is.sig & doctype=="EO"),
							numsigPR = sum(is.sig & doctype=="PR"),
							numsigMM = sum(is.sig & doctype=="MM"))
                              
significant_annual_type = unilateral_preds %>%
					mutate(is.sig = Significance_rf >= 0.355) %>%
			 		group_by(year, doctype) %>% summarize(numsig = sum(is.sig))

## distinguish unilateral activity by issue area
issues = unilateral_preds %>% 
  group_by(issue, year) %>% 
  summarize(numsig = sum(Significance_rf >= .355),
			numdocs = n()) %>%
  mutate(numsig = replace_na(numsig, 0),
		numdocs = replace_na(numdocs, 0))


## Read in other data necessary for analysis
noncere_EOs_1905_2013 = read.dta13("noncere_EOs_1905_2013.dta")
sig_EOs = read.dta13("sig EOs.dta") 
memos = read.dta13("memos.dta")
recent_noncere_EOs = read.csv("recent_noncere_EOs.csv")
mip = read.csv("US-Public-Gallups_Most_Important_Problem-21.2.csv")


## prepare "most important problem" data for merging
mip$issue[mip$majortopic==1] = "Economic"
mip$issue[mip$majortopic==2] = "Civil rights"
mip$issue[mip$majortopic==3] = "Health"
mip$issue[mip$majortopic==4] = "Agriculture"
mip$issue[mip$majortopic==5] = "Labor"
mip$issue[mip$majortopic==6] = "Education"
mip$issue[mip$majortopic==7] = "Environment"
mip$issue[mip$majortopic==8] = "Energy"
mip$issue[mip$majortopic==9] = "Immigration"
mip$issue[mip$majortopic==10] = "Transportation"
mip$issue[mip$majortopic==12] = "Law/Crime"
mip$issue[mip$majortopic==13] = "Social welfare"
mip$issue[mip$majortopic==14] = "Housing"
mip$issue[mip$majortopic==15] = "Domestic commerce"
mip$issue[mip$majortopic==16] = "Defense"
mip$issue[mip$majortopic==17] = "Technology"
mip$issue[mip$majortopic==18] = "Trade"
mip$issue[mip$majortopic==19] = "International affairs"
mip$issue[mip$majortopic==20] = "Government operations"
mip$issue[mip$majortopic==21] = "Public lands"

mip = mip %>% filter(majortopic!=25)

mip_issues <- merge(issues, mip, 
              by = c("year","issue"),all.x=TRUE,all.y=TRUE)

mip_issues = mip_issues %>% 
  filter(year >= 1956 & year<2021) %>%
  mutate(numsig = ifelse(is.na(numsig), 0, numsig)) %>%
  mutate(numdocs = ifelse(is.na(numdocs), 0, numdocs))


# Set colors for plotting in grayscale
colors3 = c("#CCCCCC","#959595", "#000000")


################################################
## Figure 1a: Number of significant documents ##
################################################

fig1a = ggplot(significant_sum, aes(x = year, y = numsig)) +
  		geom_line() +
		theme_bw(base_size=13) + ylab("Number of significant directives") + xlab("") +
		scale_x_continuous(breaks = c(1880,1900,1920,1940,1960,1980,2000,2020)) +
		theme(axis.ticks.x=element_blank()) +
		theme(axis.ticks.y=element_blank()) +
		theme(axis.title.y = element_text(margin = margin(t = 0, r = 10, b = 0, l = 0))) +
		theme(panel.grid.minor = element_blank(),panel.grid.major = element_blank())

ggsave("fig1a.pdf", fig1a, height = 5, width = 8)


###############################################
## Figure 1b: Comparison with other measures ##
###############################################

fig1b = ggplot(noncere_EOs_1905_2013, aes(x = year, y = allnoncerm_eo)) +
	  	geom_line(lwd=1,color="DARKGRAY") +
		theme_bw(base_size=13) + ylab("Number of directives") + xlab("") +
		scale_x_continuous(breaks = c(1880,1900,1920,1940,1960,1980,2000,2020)) +
		theme(axis.ticks.x=element_blank()) +
		theme(axis.ticks.y=element_blank()) +
 		geom_line(data=recent_noncere_EOs,aes(x=year,y=numdocs),color="DARKGRAY",lwd=1) +
		geom_line(data=sig_EOs,aes(x=year,y=sig_eo),lwd=.5,linetype="dotted") +
		geom_line(data=significant_sum,aes(x=year, y=numsig),lwd=1) + 
		annotate("text", x = 2009, y = 210, size = 3,label = "Significant \n directives") +
		annotate("text", x = 1950, y = 450, size = 3,label = "Nonceremonial \n executive orders \n (Bolton and Thrower)",color="gray45") +
		annotate("text", x = 1924.5, y = 5, size = 2.75,label = "Significant executive orders (Howell)") +
		theme(axis.title.y = element_text(margin = margin(t = 0, r = 10, b = 0, l = 0))) +
		theme(panel.grid.minor = element_blank(),panel.grid.major = element_blank())

ggsave("fig1b.pdf", fig1b, height = 5, width = 8)

#################################################
## Figure 1c: Variation across directive types ##
#################################################

fig1c = ggplot(significant_annual_type, aes(x=year, group=doctype, y=numsig)) + geom_line(aes(col=doctype), linewidth=1) +
		ylab("Number of significant directives") + xlab("") + theme_minimal(base_size=13) + 
		theme(legend.title=element_blank()) + scale_x_continuous(breaks = c(1880,1900,1920,1940,1960,1980,2000,2020)) +
		theme_bw(base_size=13) +
		theme(legend.position="bottom") +
		scale_color_manual(name="",values=colors3,
			breaks=c("EO", "MM","PR"),
			labels=c("Executive orders          ","Memoranda            ","Proclamations     ")) +
	theme(legend.title=element_blank()) +
	theme(axis.ticks.x=element_blank()) +
	theme(axis.ticks.y=element_blank()) +
	theme(axis.title.y = element_text(margin = margin(t = 0, r = 10, b = 0, l = 0))) +
	theme(panel.grid.minor = element_blank(),panel.grid.major = element_blank())

ggsave("fig1c.pdf", fig1c, height = 5, width = 8)

##################################################
## Figure 2: Significant activity across issues ##
##################################################

fig2 = unilateral_preds %>% 
	  group_by(issue, year) %>% 
	  summarize(numsig = sum(Significance_rf >= .355)) %>%
	  right_join(expand_grid(issue = unique(unilateral_preds$issue),
                         year= 1877:2020)) %>%
	  mutate(numsig = replace_na(numsig, 0)) %>%
	  ggplot(aes(x=year, y=numsig, group = issue)) + 
	  geom_line() +
	  facet_wrap(~issue, scales = "free_y",ncol=3) +
	  theme(strip.text = element_text(size = 7), panel.spacing = unit(.5, 'pt'))  + theme_bw(base_size=9) +
	  labs(title = "", x = "Year", y = "Number of significant directives") +
	  theme(strip.text = element_text(size = 7)) +
	  theme_bw() +
	  theme(plot.title = element_text(hjust = 0.5)) +
	  scale_y_continuous(breaks = scales::pretty_breaks(3), limits = c(0, NA)) +
	  theme(panel.grid.minor = element_blank(),panel.grid.major = element_blank())

ggsave("fig2.pdf", fig2, height = 8, width = 7)

#########################################################
## Figure 3: Directive substitution across issue areas ##
#########################################################

par(mfrow=c(1,3), mar=c(3,4,3,2))
fit_eo=lm(doctype=="EO" ~ poly(year,3,raw=TRUE), significant_preds)
fit_pr=lm(doctype=="PR" ~ poly(year,3,raw=TRUE), significant_preds)
fit_mm=lm(doctype=="MM" ~ poly(year,3,raw=TRUE), significant_preds)
plot(NA,xlim=c(1877, 2020), ylim=c(0,1), xlab="Year",
     ylab="Proportion of directives",
     main="All policy areas")
lines(1877:2020, predict(fit_eo,data.frame(year=1877:2020)), lwd=2, lty=1)
lines(1877:2020, predict(fit_pr, data.frame(year=1877:2020)), lwd=2, lty=2)
lines(1877:2020, predict(fit_mm, data.frame(year=1877:2020)), lwd=2, lty=3)
legend("topleft", c("Executive orders", "Proclamations", "Memoranda"), lty = 1:3,bty='n')

fit_eo=lm(doctype=="EO" ~ poly(year,3,raw=TRUE), significant_preds[significant_preds$policyarea_1 %in% c(1,2,3,4,5,6,7,8,10,11,12,13,14,15,21),])
fit_pr=lm(doctype=="PR" ~ poly(year,3,raw=TRUE), significant_preds[significant_preds$policyarea_1 %in% c(1,2,3,4,5,6,7,8,10,11,12,13,14,15,21),])
fit_mm=lm(doctype=="MM" ~ poly(year,3,raw=TRUE), significant_preds[significant_preds$policyarea_1 %in% c(1,2,3,4,5,6,7,8,10,11,12,13,14,15,21),])
plot(NA,xlim=c(1877, 2020), ylim=c(0,1), xlab="Year",
     ylab="Proportion of directives",
     main="Domestic Policy")
lines(1877:2020, predict(fit_eo, data.frame(year=1877:2020)), lwd=2, lty=1)
lines(1877:2020, predict(fit_pr, data.frame(year=1877:2020)), lwd=2, lty=2)
lines(1877:2020, predict(fit_mm, data.frame(year=1877:2020)), lwd=2, lty=3)
legend("topleft", c("Executive orders", "Proclamations", "Memoranda"), lty = 1:3,bty='n')

fit_eo=lm(doctype=="EO" ~ poly(year,3,raw=TRUE), significant_preds[significant_preds$policyarea_1 %in% c(9,16,18,19),])
fit_pr=lm(doctype=="PR" ~ poly(year,3,raw=TRUE), significant_preds[significant_preds$policyarea_1 %in% c(9,16,18,19),])
fit_mm=lm(doctype=="MM" ~ poly(year,3,raw=TRUE), significant_preds[significant_preds$policyarea_1 %in% c(9,16,18,19),])
plot(NA,xlim=c(1877, 2020), ylim=c(0,1), xlab="Year",
     ylab="Proportion of directives",
     main="Foreign Policy")
lines(1877:2020, predict(fit_eo, data.frame(year=1877:2020)), lwd=2, lty=1)
lines(1877:2020, predict(fit_pr, data.frame(year=1877:2020)), lwd=2, lty=2)
lines(1877:2020, predict(fit_mm, data.frame(year=1877:2020)), lwd=2, lty=3)
legend("topleft", c("Executive orders", "Proclamations", "Memoranda"), lty = 1:3,bty='n')
dev.copy(pdf,"fig3.pdf",height=5,width=13)
dev.off()

#################################
## Figure 4: Structural change ##
#################################

# First, identify changepoints

directives <- ts(significant_sum %>% pull(numsig), start=c(1877), end=c(2020), frequency=1) 
bp.directives <- breakpoints(directives ~ 1, h=4)
summary(bp.directives)

# Identify confidence intervals

ci.bp <- confint(bp.directives,het.reg=T, het.err=F)
ci.bp


# Then, plot annual number of significant directives, use vertical lines to denote
# changepoints, and use horizontal lines to indicate the mean number of annual directives
# within each period

fig4 = ggplot(significant_sum, aes(x=year, y=numsig)) + 
		geom_point(color="grey65") + theme_bw() +
		ylab("Number of Significant directives") + xlab("") + theme_minimal(base_size=13) + 
		theme(legend.title=element_blank()) + scale_x_continuous(breaks = c(1880,1900,1920,1940,1960,1980,2000,2020)) +
		theme_bw(base_size=13) +
		theme(legend.position="bottom") +
		theme(legend.title=element_blank()) +
		theme(axis.ticks.x=element_blank()) +
		theme(axis.ticks.y=element_blank()) +
		theme(axis.title.y = element_text(margin = margin(t = 0, r = 10, b = 0, l = 0))) +
		theme(panel.grid.minor = element_blank(),panel.grid.major = element_blank()) +
		geom_vline(xintercept=1903.5, linetype="dashed") +
		geom_vline(xintercept=1908.5, linetype="dashed") +
		geom_vline(xintercept=1917.5, linetype="dashed") +
		geom_vline(xintercept=1990.5, linetype="dashed") +
		geom_vline(xintercept=2015.5, linetype="dashed") +
		geom_segment(aes(x = 1877, y = mean(significant_sum$numsig[significant_sum$year>=1877 & significant_sum$year<1904]), 
			xend = 1903.5, yend = mean(significant_sum$numsig[significant_sum$year>=1877 & significant_sum$year<1904])),lwd=1.5) +
		geom_segment(aes(x = 1903.5, y = mean(significant_sum$numsig[significant_sum$year>=1904 & significant_sum$year<1909]),
			xend = 1908.5, yend = mean(significant_sum$numsig[significant_sum$year>=1904 & significant_sum$year<1909])),lwd=1.5) +
		geom_segment(aes(x = 1908.5, y = mean(significant_sum$numsig[significant_sum$year>=1909 & significant_sum$year<1918]),
			xend = 1917.5, yend = mean(significant_sum$numsig[significant_sum$year>=1909 & significant_sum$year<1918])),lwd=1.5) +
		geom_segment(aes(x = 1917.5, y = mean(significant_sum$numsig[significant_sum$year>=1918 & significant_sum$year<1991]), 
			xend = 1990.5, yend = mean(significant_sum$numsig[significant_sum$year>=1918 & significant_sum$year<1991])),lwd=1.5) +
		geom_segment(aes(x = 1990.5, y = mean(significant_sum$numsig[significant_sum$year>=1991 & significant_sum$year<2016]), 
			xend = 2015.5, yend = mean(significant_sum$numsig[significant_sum$year>=1991 & significant_sum$year<2016])),lwd=1.5) +
		geom_segment(aes(x = 2015.5, y = mean(significant_sum$numsig[significant_sum$year>=2016]), 
			xend = 2021, yend = mean(significant_sum$numsig[significant_sum$year>=2016])),lwd=1.5)
ggsave("fig4.pdf", fig4, height = 5, width = 8)


########################################################
## Figure 5:					      ##
## Correlation between unilateral action and MIP data ##
## (replication of Jones et al 2009)		      ##
########################################################


# First, estimate correlations

## Only significant documents ##
cor.test(mip_issues$percent[mip_issues$issue=="Economic"], mip_issues$numsig[mip_issues$issue=="Economic"], use="complete.obs")
cor.test(mip_issues$percent[mip_issues$issue=="Civil rights"], mip_issues$numsig[mip_issues$issue=="Civil rights"], use="complete.obs")
cor.test(mip_issues$percent[mip_issues$issue=="Health"], mip_issues$numsig[mip_issues$issue=="Health"], use="complete.obs")
cor.test(mip_issues$percent[mip_issues$issue=="Agriculture"], mip_issues$numsig[mip_issues$issue=="Agriculture"], use="complete.obs")
cor.test(mip_issues$percent[mip_issues$issue=="Labor"], mip_issues$numsig[mip_issues$issue=="Labor"], use="complete.obs")
cor.test(mip_issues$percent[mip_issues$issue=="Education"], mip_issues$numsig[mip_issues$issue=="Education"], use="complete.obs")
cor.test(mip_issues$percent[mip_issues$issue=="Environment"], mip_issues$numsig[mip_issues$issue=="Environment"], use="complete.obs")
cor.test(mip_issues$percent[mip_issues$issue=="Energy"], mip_issues$numsig[mip_issues$issue=="Energy"], use="complete.obs")
cor.test(mip_issues$percent[mip_issues$issue=="Immigration"], mip_issues$numsig[mip_issues$issue=="Immigration"], use="complete.obs")
cor.test(mip_issues$percent[mip_issues$issue=="Transportation"], mip_issues$numsig[mip_issues$issue=="Transportation"], use="complete.obs")
cor.test(mip_issues$percent[mip_issues$issue=="Law/Crime"], mip_issues$numsig[mip_issues$issue=="Law/Crime"], use="complete.obs")
cor.test(mip_issues$percent[mip_issues$issue=="Social welfare"], mip_issues$numsig[mip_issues$issue=="Social welfare"], use="complete.obs")
cor.test(mip_issues$percent[mip_issues$issue=="Housing"], mip_issues$numsig[mip_issues$issue=="Housing"], use="complete.obs")
cor.test(mip_issues$percent[mip_issues$issue=="Domestic commerce"], mip_issues$numsig[mip_issues$issue=="Domestic commerce"], use="complete.obs")
cor.test(mip_issues$percent[mip_issues$issue=="Defense"], mip_issues$numsig[mip_issues$issue=="Defense"], use="complete.obs")
cor.test(mip_issues$percent[mip_issues$issue=="Technology"], mip_issues$numsig[mip_issues$issue=="Technology"], use="complete.obs")
cor.test(mip_issues$percent[mip_issues$issue=="Trade"], mip_issues$numsig[mip_issues$issue=="Trade"], use="complete.obs")
cor.test(mip_issues$percent[mip_issues$issue=="International affairs"], mip_issues$numsig[mip_issues$issue=="International affairs"], use="complete.obs")
cor.test(mip_issues$percent[mip_issues$issue=="Government operations"], mip_issues$numsig[mip_issues$issue=="Government operations"], use="complete.obs")
cor.test(mip_issues$percent[mip_issues$issue=="Public lands"], mip_issues$numsig[mip_issues$issue=="Public lands"], use="complete.obs")


## All documents
cor.test(mip_issues$percent[mip_issues$issue=="Economic"], mip_issues$numdocs[mip_issues$issue=="Economic"], use="complete.obs")
cor.test(mip_issues$percent[mip_issues$issue=="Civil rights"], mip_issues$numdocs[mip_issues$issue=="Civil rights"], use="complete.obs")
cor.test(mip_issues$percent[mip_issues$issue=="Health"], mip_issues$numdocs[mip_issues$issue=="Health"], use="complete.obs")
cor.test(mip_issues$percent[mip_issues$issue=="Agriculture"], mip_issues$numdocs[mip_issues$issue=="Agriculture"], use="complete.obs")
cor.test(mip_issues$percent[mip_issues$issue=="Labor"], mip_issues$numdocs[mip_issues$issue=="Labor"], use="complete.obs")
cor.test(mip_issues$percent[mip_issues$issue=="Education"], mip_issues$numdocs[mip_issues$issue=="Education"], use="complete.obs")
cor.test(mip_issues$percent[mip_issues$issue=="Environment"], mip_issues$numdocs[mip_issues$issue=="Environment"], use="complete.obs")
cor.test(mip_issues$percent[mip_issues$issue=="Energy"], mip_issues$numdocs[mip_issues$issue=="Energy"], use="complete.obs")
cor.test(mip_issues$percent[mip_issues$issue=="Immigration"], mip_issues$numdocs[mip_issues$issue=="Immigration"], use="complete.obs")
cor.test(mip_issues$percent[mip_issues$issue=="Transportation"], mip_issues$numdocs[mip_issues$issue=="Transportation"], use="complete.obs")
cor.test(mip_issues$percent[mip_issues$issue=="Law/Crime"], mip_issues$numdocs[mip_issues$issue=="Law/Crime"], use="complete.obs")
cor.test(mip_issues$percent[mip_issues$issue=="Social welfare"], mip_issues$numdocs[mip_issues$issue=="Social welfare"], use="complete.obs")
cor.test(mip_issues$percent[mip_issues$issue=="Housing"], mip_issues$numdocs[mip_issues$issue=="Housing"], use="complete.obs")
cor.test(mip_issues$percent[mip_issues$issue=="Domestic commerce"], mip_issues$numdocs[mip_issues$issue=="Domestic commerce"], use="complete.obs")
cor.test(mip_issues$percent[mip_issues$issue=="Defense"], mip_issues$numdocs[mip_issues$issue=="Defense"], use="complete.obs")
cor.test(mip_issues$percent[mip_issues$issue=="Technology"], mip_issues$numdocs[mip_issues$issue=="Technology"], use="complete.obs")
cor.test(mip_issues$percent[mip_issues$issue=="Trade"], mip_issues$numdocs[mip_issues$issue=="Trade"], use="complete.obs")
cor.test(mip_issues$percent[mip_issues$issue=="International affairs"], mip_issues$numdocs[mip_issues$issue=="International affairs"], use="complete.obs")
cor.test(mip_issues$percent[mip_issues$issue=="Government operations"], mip_issues$numdocs[mip_issues$issue=="Government operations"], use="complete.obs")
cor.test(mip_issues$percent[mip_issues$issue=="Public lands"], mip_issues$numdocs[mip_issues$issue=="Public lands"], use="complete.obs")


# Then, plot correlations in descending order by issue area

sig = c(.41,.4,.39,.36,.33,.33,.32,.31,.28,.28,.25,.17,.11,.05,.03,.01,-.08,-.16,-.16,-.29)
sig_low = c(.19,.18,.16,.12,.1,.09,.08,.07,.04,.04,.01,-.07,-.13,-.19,-.21,-.23,-.31,-.39,-.39,-.49)
sig_hi = c(.6,.59,.58,.55,.53,.53,.52,.51,.49,.49,.47,.4,.35,.29,.28,.26,.17,.09,.08,-.04)
sig_issue = c("Health","Immigration","Civil rights","Agriculture","Tranportation",
			"Govt operations","Education","Labor","Technology","Trade","Energy",
			"Law/Crime","Domestic commerce","Housing","Environment","Economic",
			"Social welfare","Public lands","International affairs","Defense")
mean(sig)

all = c(.58,.55,.28,.27,.24,.23,.19,.15,.14,.11,.1,.06,.05,.02,.01,-.04,-.19,-.21,-.21,-.24)
all_low = c(.39,.35,.04,.02,-.01,-.01,-.06,-.1,-.11,-.14,-.15,-.19,-.2,-.23,-.23,-.28,-.42,-.43,-.45,-.45)
all_hi = c(.72,.7,.49,.48,.45,.45,.41,.38,.37,.35,.33,.3,.29,.26,.26,.21,.05,.04,.04,.01)
all_issue = c("Education","Health","Social welfare","Labor","Law/crime","Civil rights",
			"Govt operations","Immigration","Transportation","Energy","Trade",
			"Economic","Housing","Technology","Agriculture","Public lands",
			"Domestic commerce","Defense",
			"Environment","International affairs")
order = seq(from=20,to=1,by=-1)
mean(all)

jones = c(.66,.38,.29,.23,.18,.16,.16,.12,.11,.09,.08,.08,
		.03,.01,.00,.00,-.03,-.06)
jones_issues = c("Energy","Environment","Agriculture","Housing","Law/crime",
		"Health","Education","Technology","Welfare","Trade",
		"Civil rights","Economics","Govt operations",
		"Public lands","Defense","Commerce","Labor","Transportation")	
jones_order = seq(from=18,to=1,by=-1)

par(mar=c(3.5,7,3,1), mgp=c(2,.75,0),mfrow=c(1,3))
plot(jones,jones_order,xlim=c(-.7,.8),ylim=c(0.75,18.25),cex.lab=1,cex.axis=1,pch=20,cex=2,ylab="",
	xaxt="n",yaxt="n",xlab="Correlation coefficient",main="Executive orders, 1956-2002\n (from Jones et al. 2009)",
	col=c("BLACK","BLACK","BLACK","GREY","GREY","GREY","GREY","GREY","GREY","GREY","GREY",
		"GREY","GREY","GREY","GREY","GREY","GREY","GREY"))
abline(v=0,lty=2)
par(mgp=c(2,.1,0))
axis(2,at=jones_order,las=2,cex.axis=1,labels=jones_issues,tck=0)
axis(1,tck=0)

par(mar=c(3.5,7,3,1), mgp=c(2,.75,0))
plot(all,order,xlim=c(-.7,.8),ylim=c(0.75,20.25),cex.lab=1,cex.axis=1,pch=20,cex=2,ylab="",
	xaxt="n",yaxt="n",xlab="Correlation coefficient",main="All directives, 1956-2020",
		col=c("BLACK","BLACK","BLACK","BLACK","GREY","GREY","GREY","GREY","GREY","GREY","GREY",
		"GREY","GREY","GREY","GREY","GREY","GREY","GREY","GREY","GREY"))
abline(v=0,lty=2)
segments(all_low,order,all_hi,order,col=c("BLACK","BLACK","BLACK","BLACK","GREY","GREY","GREY","GREY","GREY","GREY","GREY",
		"GREY","GREY","GREY","GREY","GREY","GREY","GREY","GREY","GREY"))
par(mgp=c(2,.1,0))
axis(2,at=order,las=2,cex.axis=1,labels=all_issue,tck=0)
axis(1,tck=0)

par(mar=c(3.5,7,3,1), mgp=c(2,.75,0))
plot(sig,order,xlim=c(-.7,.8),ylim=c(0.75,20.25),cex.lab=1,cex.axis=1,pch=20,cex=2,ylab="",
	xaxt="n",yaxt="n",xlab="Correlation coefficient",main="Significant directives, 1956-2020",
	col=c("BLACK","BLACK","BLACK","BLACK","BLACK","BLACK","BLACK","BLACK","BLACK","BLACK","BLACK","GREY","GREY",
		"GREY","GREY","GREY","GREY","GREY","GREY","BLACK"))
abline(v=0,lty=2)
segments(sig_low,order,sig_hi,order,col=c("BLACK","BLACK","BLACK","BLACK","BLACK","BLACK","BLACK","BLACK","BLACK","BLACK","BLACK","GREY","GREY",
		"GREY","GREY","GREY","GREY","GREY","GREY","BLACK"))
par(mgp=c(2,.1,0))
axis(2,at=order,las=2,cex.axis=1,labels=sig_issue,tck=0)
axis(1,tck=0)

dev.copy(pdf,"fig5.pdf",height=6,width=12)
dev.off()



#################################################
## Figure A.1: All directives by year and type ##
#################################################

figa1 = ggplot(unilateral_annual_type, aes(x=year, group=doctype, y=numdocs)) + geom_line(aes(col=doctype), size=1) +
		ylab("Number of Directives") + xlab("") + theme_minimal(base_size=13) + 
		theme(legend.title=element_blank()) + scale_x_continuous(breaks = c(1880,1900,1920,1940,1960,1980,2000,2020)) +
		theme_bw(base_size=13) +
		theme(legend.position="bottom") +
		scale_color_manual(name="",values=colors3,
                         breaks=c("EO", "MM","PR"),
                         labels=c("Executive orders          ","Memoranda            ","Proclamations     ")) +
		theme(legend.title=element_blank()) +
		theme(axis.ticks.x=element_blank()) +
		theme(axis.ticks.y=element_blank()) +
		theme(axis.title.y = element_text(margin = margin(t = 0, r = 10, b = 0, l = 0))) +
		theme(panel.grid.minor = element_blank(),panel.grid.major = element_blank())
ggsave("figa1.pdf", figa1, height = 5, width = 8)

########################################################
## Figure A.2: Distribution of significance estimates ##
########################################################

figa2 = ggplot(unilateral_preds, aes(x=Significance_rf)) + geom_density(fill="LIGHTGRAY") + ylab("Density") +
  		theme_bw(base_size=13) + theme(legend.title=element_blank()) + 
		xlab("Significance estimates") +
 		theme(axis.title.y = element_text(margin = margin(t = 0, r = 10, b = 0, l = 0))) +
		theme(panel.grid.minor = element_blank(),panel.grid.major = element_blank())
ggsave("figa2.pdf", figa2, height = 5, width = 8)


###########################################
## Figure B.1: Temporal Cross-validation ##
###########################################

load("train_test_data.RData")
set.seed(10012)

idx = sample(x=1:10, size=nrow(chiou2), replace=T)
cross_val_rf = function(x, thresh){
  ytest = as.numeric(chiou2$Significance>thresh)
  train = chiou_subset[idx!=x,]
  train = quanteda::convert(train, to="data.frame")
  test = chiou_subset[idx==x,]
  test = quanteda::convert(test, to="data.frame")
  train$outcome = factor(ytest[idx!=x])
  m = ranger::ranger(dependent.variable.name="outcome", data=train,
                     probability=FALSE, write.forest=TRUE, num.trees=2000)
  preds = predict(m, test, predict.all=TRUE)
  p2 = rowMeans(preds$predictions)
  out = cor(p2, ytest[idx == x]) 
  print(out)
  out2 = list(predvec = p2, truth = ytest[idx==x])
  return(out2)
}

get_roc = function(thresh=.5){
  cv_ests = sapply(1:10, FUN=function(x) cross_val_rf(x, thresh)) 
  preds = unlist(cv_ests[c(1,3,5,7,9,11,13,15,17,19)])
  truth =  unlist(cv_ests[c(2,4,6,8,10,12,14,16,18,20)])
  out = roc(response=truth, predictor = preds)
  return(out)
}
cv_roc = get_roc()
cv_roc # 0.9056

idx = cut(1:nrow(chiou_subset), 10)
cv_ests = sapply(unique(idx), FUN=function(x) cross_val_rf(x, 0.5)) 
preds = unlist(cv_ests[c(1,3,5,7,9,11,13,15,17)])
truth =  unlist(cv_ests[c(2,4,6,8,10,12,14,16,18)])
temporal_cv_roc = roc(response = truth, predictor=preds)
temporal_cv_roc # .884


png(filename = "temporal_auc_plots.png")
plot(cv_roc)
plot(temporal_cv_roc, lty=2, add=T)
# Then add text and line segments as legend
segments(x0=0.5, x1=0.4, y0=0.2, y1=0.2, lw=2)
segments(x0=0.5, x1=0.4, y0=0.125, y1=0.125, lty=2, lwd=2)
text(x=0.115, y=.2, "Normal AUC = 0.906")
text(x=0.135, y=.125, "Temporal AUC = 0.884")
dev.off()


############################################
## Table B.1				  ##
## Feature importance, significance model ##
############################################

load("train_test_data.RData")
set.seed(10012)

temp = quanteda::convert(chiou_subset, to="data.frame")
temp$outcome=factor(as.numeric(chiou2$Significance>.5))

dim(tdm1) #101186, 30229
tdm1.1 = quanteda::convert(tdm1[1:50000,], to="matrix")
tdm1.2 =  quanteda::convert(tdm1[50001:nrow(tdm1),], to="matrix")

tdm1 = rbind(tdm1.1,tdm1.2)

rm(tdm1.1)
rm(tdm1.2)
gc()

m = ranger(dependent.variable.name="outcome", data=temp[,-1],
           probability=TRUE, write.forest=FALSE, num.trees=5000,
           importance = "impurity_corrected")

imp = importance(m)
imp = imp[order(imp, decreasing=T)]
imp = round(imp, 3)
imp3 = cbind(names(imp), imp)


imp3 = imp3[nchar(imp3[,1])>6,]
imp3 = imp3[!grepl("_", imp3[,1]),]


print(xtable(imp3[1:30,]), include.rownames=FALSE)

tab2 = cbind(imp3[1:30,], imp3[2701:2730,])
print(xtable(tab2, digits=3), include.rownames=FALSE,
      file="feature_importance_significance.tex")


#####################################
## Table B.2			   ##
## Feature importance, policy area ##
#####################################

load("tmp_policy_classifier.RData")

set.seed(10012)
m = ranger::ranger(dependent.variable.name="train.y2", data=train.x2,
                   probability=FALSE, write.forest=FALSE, num.trees=5000,
                   importance = "impurity_corrected")

imp = importance(m)
imp = imp[order(imp, decreasing=T)]
imp = round(imp, 3)
imp3 = cbind(names(imp), imp)
imp3 = imp3[nchar(imp3[,1])>6,]
imp3 = imp3[!grepl("_", imp3[,1]),]

print(xtable(imp3[1:30,]), include.rownames=FALSE)

tab2 = cbind(imp3[1:30,], imp3[c(2601:2630),])
print(xtable(tab2, digits=3), include.rownames=FALSE,
      file="feature_importance_policy.tex")


########################################
## Table B.3			      ##
## Descriptive stats on proclamations ##
########################################

# Proclamation numbers and titles from 
# https://www.federalregister.gov/presidential-documents/proclamations/donald-trump/2019
# (accessed 13 October 2023)
# Then, add significance estimates for each 

round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9836" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9837" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9838" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9839" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9840" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9841" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9842" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9843" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9844" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9845" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9846" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9847" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9848" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9849" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9850" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9851" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9852" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9853" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9854" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9855" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9856" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9857" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9858" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9859" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9860" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9861" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9862" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9863" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9864" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9865" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9866" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9867" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9868" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9869" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9870" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9871" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9872" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9873" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9874" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9875" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9876" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9877" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9878" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9879" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9880" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9881" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9882" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9883" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9884" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9885" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9886" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9887" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9888" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9889" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9890" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9891" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9892" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9893" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9894" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9895" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9896" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9897" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9898" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9899" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9900" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9901" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9902" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9903" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9904" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9905" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9906" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9907" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9908" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9909" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9910" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9911" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9912" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9913" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9914" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9915" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9916" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9917" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9918" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9919" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9920" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9921" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9922" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9923" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9924" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9925" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9926" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9927" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9928" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9929" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9930" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9931" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9932" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9933" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9934" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9935" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9936" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9937" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9938" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9939" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9940" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9941" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9942" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9943" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9944" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9945" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9946" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9947" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9948" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9949" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9950" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9951" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9952" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9953" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9954" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9955" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9956" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9957" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9958" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9959" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9960" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9961" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9962" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9963" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9964" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9965" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9966" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9967" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9968" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9969" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9970" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9971" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9972" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9973" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9974" & unilateral_preds$traintest=="test"],digits=2)
round(unilateral_preds$Significance_rf[unilateral_preds$accession=="2019-PR-9975" & unilateral_preds$traintest=="test"],digits=2)

# Then, add estimated issue area

unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9836" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9837" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9838" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9839" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9840" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9841" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9842" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9843" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9844" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9845" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9846" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9847" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9848" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9849" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9850" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9851" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9852" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9853" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9854" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9855" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9856" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9857" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9858" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9859" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9860" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9861" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9862" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9863" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9864" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9865" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9866" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9867" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9868" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9869" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9870" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9871" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9872" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9873" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9874" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9875" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9876" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9877" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9878" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9879" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9880" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9881" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9882" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9883" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9884" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9885" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9886" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9887" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9888" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9889" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9890" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9891" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9892" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9893" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9894" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9895" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9896" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9897" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9898" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9899" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9900" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9901" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9902" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9903" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9904" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9905" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9906" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9907" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9908" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9909" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9910" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9911" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9912" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9913" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9914" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9915" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9916" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9917" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9918" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9919" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9920" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9921" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9922" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9923" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9924" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9925" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9926" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9927" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9928" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9929" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9930" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9931" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9932" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9933" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9934" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9935" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9936" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9937" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9938" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9939" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9940" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9941" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9942" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9943" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9944" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9945" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9946" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9947" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9948" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9949" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9950" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9951" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9952" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9953" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9954" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9955" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9956" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9957" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9958" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9959" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9960" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9961" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9962" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9963" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9964" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9965" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9966" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9967" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9968" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9969" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9970" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9971" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9972" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9973" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9974" & unilateral_preds$traintest=="test"]
unilateral_preds$issue[unilateral_preds$accession=="2019-PR-9975" & unilateral_preds$traintest=="test"]

# Now create an indicator to distinguish the 15 nonceremonial proclamations from the others

unilateral_preds$ceremonial = ifelse(
    ((unilateral_preds$accession %in% c("2019-PR-9842","2019-PR-9844",
						"2019-PR-9852","2019-PR-9880",
						"2019-PR-9886","2019-PR-9887",
						"2019-PR-9888","2019-PR-9893",
						"2019-PR-9894","2019-PR-9902",
						"2019-PR-9931","2019-PR-9932",
						"2019-PR-9845","2019-PR-9855",
						"2019-PR-9874")) &
        (unilateral_preds$year == 2019) &
        (unilateral_preds$doctype == "PR") &
	  (unilateral_preds$traintest=="test")),
    0,  
    1)

unilateral_preds$ceremonial = ifelse(
    ((unilateral_preds$year == 2019) &
        (unilateral_preds$doctype == "PR") &
	  (unilateral_preds$traintest=="test")),
    unilateral_preds$ceremonial,
    NA)

# Calculate the mean significance for ceremonial and nonceremonial proclamations
table(unilateral_preds$ceremonial)
summary(unilateral_preds$Significance_rf[unilateral_preds$ceremonial==0])
summary(unilateral_preds$Significance_rf[unilateral_preds$ceremonial==1])


########################################################
## Appendix C: Other linear/polynomial specifications ##
########################################################

# Figure C.1a

par(mfrow=c(1,3), mar=c(3,4,3,2))
fit_eo=lm(doctype=="EO" ~ poly(year,1,raw=TRUE), significant_preds)
fit_pr=lm(doctype=="PR" ~ poly(year,1,raw=TRUE), significant_preds)
fit_mm=lm(doctype=="MM" ~ poly(year,1,raw=TRUE), significant_preds)
plot(NA,xlim=c(1877, 2020), ylim=c(0,1), xlab="Year",
     ylab="Proportion of directives",
     main="All policy areas")
lines(1877:2020, predict(fit_eo, data.frame(year=1877:2020)), lwd=2, lty=1)
lines(1877:2020, predict(fit_pr, data.frame(year=1877:2020)), lwd=2, lty=2)
lines(1877:2020, predict(fit_mm, data.frame(year=1877:2020)), lwd=2, lty=3)
legend("topleft", c("Executive orders", "Proclamations", "Memoranda"), lty = 1:3,bty='n')

fit_eo=lm(doctype=="EO" ~ poly(year,1,raw=TRUE), significant_preds[significant_preds$policyarea_1 %in% c(1,2,3,4,5,6,7,8,10,11,12,13,14,15,21),])
fit_pr=lm(doctype=="PR" ~ poly(year,1,raw=TRUE), significant_preds[significant_preds$policyarea_1 %in% c(1,2,3,4,5,6,7,8,10,11,12,13,14,15,21),])
fit_mm=lm(doctype=="MM" ~ poly(year,1,raw=TRUE), significant_preds[significant_preds$policyarea_1 %in% c(1,2,3,4,5,6,7,8,10,11,12,13,14,15,21),])
plot(NA,xlim=c(1877, 2020), ylim=c(0,1), xlab="Year",
     ylab="Proportion of directives",
     main="Domestic Policy")
lines(1877:2020, predict(fit_eo, data.frame(year=1877:2020)), lwd=2, lty=1)
lines(1877:2020, predict(fit_pr, data.frame(year=1877:2020)), lwd=2, lty=2)
lines(1877:2020, predict(fit_mm, data.frame(year=1877:2020)), lwd=2, lty=3)
legend("topleft", c("Executive orders", "Proclamations", "Memoranda"), lty = 1:3,bty='n')

fit_eo=lm(doctype=="EO" ~ poly(year,1,raw=TRUE), significant_preds[significant_preds$policyarea_1 %in% c(9,16,18,19),])
fit_pr=lm(doctype=="PR" ~ poly(year,1,raw=TRUE), significant_preds[significant_preds$policyarea_1 %in% c(9,16,18,19),])
fit_mm=lm(doctype=="MM" ~ poly(year,1,raw=TRUE), significant_preds[significant_preds$policyarea_1 %in% c(9,16,18,19),])
plot(NA,xlim=c(1877, 2020), ylim=c(0,1), xlab="Year",
     ylab="Proportion of directives",
     main="Foreign Policy")
lines(1877:2020, predict(fit_eo, data.frame(year=1877:2020)), lwd=2, lty=1)
lines(1877:2020, predict(fit_pr, data.frame(year=1877:2020)), lwd=2, lty=2)
lines(1877:2020, predict(fit_mm, data.frame(year=1877:2020)), lwd=2, lty=3)
legend("topleft", c("Executive orders", "Proclamations", "Memoranda"), lty = 1:3,bty='n')
dev.copy(pdf,"figca1.pdf",height=5,width=13)
dev.off()


# Figure C.1b

par(mfrow=c(1,3), mar=c(3,4,3,2))
fit_eo=lm(doctype=="EO" ~ poly(year,2,raw=TRUE), significant_preds)
fit_pr=lm(doctype=="PR" ~ poly(year,2,raw=TRUE), significant_preds)
fit_mm=lm(doctype=="MM" ~ poly(year,2,raw=TRUE), significant_preds)
plot(NA,xlim=c(1877, 2020), ylim=c(0,1), xlab="Year",
     ylab="Proportion of directives",
     main="All policy areas")
lines(1877:2020, predict(fit_eo, data.frame(year=1877:2020)), lwd=2, lty=1)
lines(1877:2020, predict(fit_pr, data.frame(year=1877:2020)), lwd=2, lty=2)
lines(1877:2020, predict(fit_mm, data.frame(year=1877:2020)), lwd=2, lty=3)
legend("topleft", c("Executive orders", "Proclamations", "Memoranda"), lty = 1:3,bty='n')

fit_eo=lm(doctype=="EO" ~ poly(year,2,raw=TRUE), significant_preds[significant_preds$policyarea_1 %in% c(1,2,3,4,5,6,7,8,10,11,12,13,14,15,21),])
fit_pr=lm(doctype=="PR" ~ poly(year,2,raw=TRUE), significant_preds[significant_preds$policyarea_1 %in% c(1,2,3,4,5,6,7,8,10,11,12,13,14,15,21),])
fit_mm=lm(doctype=="MM" ~ poly(year,2,raw=TRUE), significant_preds[significant_preds$policyarea_1 %in% c(1,2,3,4,5,6,7,8,10,11,12,13,14,15,21),])
plot(NA,xlim=c(1877, 2020), ylim=c(0,1), xlab="Year",
     ylab="Proportion of directives",
     main="Domestic Policy")
lines(1877:2020, predict(fit_eo, data.frame(year=1877:2020)), lwd=2, lty=1)
lines(1877:2020, predict(fit_pr, data.frame(year=1877:2020)), lwd=2, lty=2)
lines(1877:2020, predict(fit_mm, data.frame(year=1877:2020)), lwd=2, lty=3)
legend("topleft", c("Executive orders", "Proclamations", "Memoranda"), lty = 1:3,bty='n')

fit_eo=lm(doctype=="EO" ~ poly(year,2,raw=TRUE), significant_preds[significant_preds$policyarea_1 %in% c(9,16,18,19),])
fit_pr=lm(doctype=="PR" ~ poly(year,2,raw=TRUE), significant_preds[significant_preds$policyarea_1 %in% c(9,16,18,19),])
fit_mm=lm(doctype=="MM" ~ poly(year,2,raw=TRUE), significant_preds[significant_preds$policyarea_1 %in% c(9,16,18,19),])
plot(NA,xlim=c(1877, 2020), ylim=c(0,1), xlab="Year",
     ylab="Proportion of directives",
     main="Foreign Policy")
lines(1877:2020, predict(fit_eo, data.frame(year=1877:2020)), lwd=2, lty=1)
lines(1877:2020, predict(fit_pr, data.frame(year=1877:2020)), lwd=2, lty=2)
lines(1877:2020, predict(fit_mm, data.frame(year=1877:2020)), lwd=2, lty=3)
legend("topleft", c("Executive orders", "Proclamations", "Memoranda"), lty = 1:3,bty='n')
dev.copy(pdf,"figc1b.pdf",height=5,width=13)
dev.off()

# Figure C.2a

par(mfrow=c(1,3), mar=c(3,4,3,2))
fit_eo=lm(doctype=="EO" ~ poly(year,4,raw=TRUE), significant_preds)
fit_pr=lm(doctype=="PR" ~ poly(year,4,raw=TRUE), significant_preds)
fit_mm=lm(doctype=="MM" ~ poly(year,4,raw=TRUE), significant_preds)
plot(NA,xlim=c(1877, 2020), ylim=c(0,1), xlab="Year",
     ylab="Proportion of directives",
     main="All policy areas")
lines(1877:2020, predict(fit_eo, data.frame(year=1877:2020)), lwd=2, lty=1)
lines(1877:2020, predict(fit_pr, data.frame(year=1877:2020)), lwd=2, lty=2)
lines(1877:2020, predict(fit_mm, data.frame(year=1877:2020)), lwd=2, lty=3)
legend("topleft", c("Executive orders", "Proclamations", "Memoranda"), lty = 1:3,bty='n')

fit_eo=lm(doctype=="EO" ~ poly(year,4,raw=TRUE), significant_preds[significant_preds$policyarea_1 %in% c(1,2,3,4,5,6,7,8,10,11,12,13,14,15,21),])
fit_pr=lm(doctype=="PR" ~ poly(year,4,raw=TRUE), significant_preds[significant_preds$policyarea_1 %in% c(1,2,3,4,5,6,7,8,10,11,12,13,14,15,21),])
fit_mm=lm(doctype=="MM" ~ poly(year,4,raw=TRUE), significant_preds[significant_preds$policyarea_1 %in% c(1,2,3,4,5,6,7,8,10,11,12,13,14,15,21),])
plot(NA,xlim=c(1877, 2020), ylim=c(0,1), xlab="Year",
     ylab="Proportion of directives",
     main="Domestic Policy")
lines(1877:2020, predict(fit_eo, data.frame(year=1877:2020)), lwd=2, lty=1)
lines(1877:2020, predict(fit_pr, data.frame(year=1877:2020)), lwd=2, lty=2)
lines(1877:2020, predict(fit_mm, data.frame(year=1877:2020)), lwd=2, lty=3)
legend("topleft", c("Executive orders", "Proclamations", "Memoranda"), lty = 1:3,bty='n')

fit_eo=lm(doctype=="EO" ~ poly(year,4,raw=TRUE), significant_preds[significant_preds$policyarea_1 %in% c(9,16,18,19),])
fit_pr=lm(doctype=="PR" ~ poly(year,4,raw=TRUE), significant_preds[significant_preds$policyarea_1 %in% c(9,16,18,19),])
fit_mm=lm(doctype=="MM" ~ poly(year,4,raw=TRUE), significant_preds[significant_preds$policyarea_1 %in% c(9,16,18,19),])
plot(NA,xlim=c(1877, 2020), ylim=c(0,1), xlab="Year",
     ylab="Proportion of directives",
     main="Foreign Policy")
lines(1877:2020, predict(fit_eo, data.frame(year=1877:2020)), lwd=2, lty=1)
lines(1877:2020, predict(fit_pr, data.frame(year=1877:2020)), lwd=2, lty=2)
lines(1877:2020, predict(fit_mm, data.frame(year=1877:2020)), lwd=2, lty=3)
legend("topleft", c("Executive orders", "Proclamations", "Memoranda"), lty = 1:3,bty='n')
dev.copy(pdf,"figc2a.pdf",height=5,width=13)
dev.off()

# Figure C.2b

par(mfrow=c(1,3), mar=c(3,4,3,2))
fit_eo=lm(doctype=="EO" ~ poly(year,5,raw=TRUE), significant_preds)
fit_pr=lm(doctype=="PR" ~ poly(year,5,raw=TRUE), significant_preds)
fit_mm=lm(doctype=="MM" ~ poly(year,5,raw=TRUE), significant_preds)
plot(NA,xlim=c(1877, 2020), ylim=c(0,1), xlab="Year",
     ylab="Proportion of directives",
     main="All policy areas")
lines(1877:2020, predict(fit_eo, data.frame(year=1877:2020)), lwd=2, lty=1)
lines(1877:2020, predict(fit_pr, data.frame(year=1877:2020)), lwd=2, lty=2)
lines(1877:2020, predict(fit_mm, data.frame(year=1877:2020)), lwd=2, lty=3)
legend("topleft", c("Executive orders", "Proclamations", "Memoranda"), lty = 1:3,bty='n')

fit_eo=lm(doctype=="EO" ~ poly(year,5,raw=TRUE), significant_preds[significant_preds$policyarea_1 %in% c(1,2,3,4,5,6,7,8,10,11,12,13,14,15,21),])
fit_pr=lm(doctype=="PR" ~ poly(year,5,raw=TRUE), significant_preds[significant_preds$policyarea_1 %in% c(1,2,3,4,5,6,7,8,10,11,12,13,14,15,21),])
fit_mm=lm(doctype=="MM" ~ poly(year,5,raw=TRUE), significant_preds[significant_preds$policyarea_1 %in% c(1,2,3,4,5,6,7,8,10,11,12,13,14,15,21),])
plot(NA,xlim=c(1877, 2020), ylim=c(0,1), xlab="Year",
     ylab="Proportion of directives",
     main="Domestic Policy")
lines(1877:2020, predict(fit_eo, data.frame(year=1877:2020)), lwd=2, lty=1)
lines(1877:2020, predict(fit_pr, data.frame(year=1877:2020)), lwd=2, lty=2)
lines(1877:2020, predict(fit_mm, data.frame(year=1877:2020)), lwd=2, lty=3)
legend("topleft", c("Executive orders", "Proclamations", "Memoranda"), lty = 1:3,bty='n')

fit_eo=lm(doctype=="EO" ~ poly(year,5,raw=TRUE), significant_preds[significant_preds$policyarea_1 %in% c(9,16,18,19),])
fit_pr=lm(doctype=="PR" ~ poly(year,5,raw=TRUE), significant_preds[significant_preds$policyarea_1 %in% c(9,16,18,19),])
fit_mm=lm(doctype=="MM" ~ poly(year,5,raw=TRUE), significant_preds[significant_preds$policyarea_1 %in% c(9,16,18,19),])
plot(NA,xlim=c(1877, 2020), ylim=c(0,1), xlab="Year",
     ylab="Proportion of directives",
     main="Foreign Policy")
lines(1877:2020, predict(fit_eo, data.frame(year=1877:2020)), lwd=2, lty=1)
lines(1877:2020, predict(fit_pr, data.frame(year=1877:2020)), lwd=2, lty=2)
lines(1877:2020, predict(fit_mm, data.frame(year=1877:2020)), lwd=2, lty=3)
legend("topleft", c("Executive orders", "Proclamations", "Memoranda"), lty = 1:3,bty='n')
dev.copy(pdf,"figc2b.pdf",height=5,width=13)
dev.off()




